#######################################################

library(data.table)
library(ggplot2)
library(pheatmap)
library(dplyr)
library(optparse)
library(GSVA)

##########################################################################################
option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--gene_list_file"), type = "character"),
    make_option(c("--out_name"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 输出文件
    input_file <- "~/20231121_singleMuti/results/qc_atac/GeneExpression.MeanByCellType.magic.tsv"

    ## 基因列表
    gene_list_file <- "~/20231121_singleMuti/config/reportDiffGene/2018_cell_stem_cell.list"

    ## 输出名字
    out_name <- "2018_cell_stem_cell"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/gsva"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
gene_list_file <- opt$gene_list_file
out_name <- opt$out_name
out_path <- opt$out_path

dir.create( out_path , recursive = T )

###########################################################################################
## 读入文件
lab_magic_RNA_meanbycluster <- read.table(input_file,header=T, sep = "\t", stringsAsFactors = FALSE)
cell_research_dynamics <- fread(gene_list_file)

###########################################################################################
## 细胞顺序
grp_order2 = c("SSC",
"Differenting&Differented SPG",
"Leptotene",
"Zygotene",
"Patchytene",
"Diplotene",
"Early stage of spermatids",
"Round&ElongateS.tids",
"Sperm",
"Leydig cells",
"Myoid cells",
"Pericytes",
"Sertoli cells",
"Endothelial cells",
"NKT cells",
"Macrophages")

grp_order2 <- gsub( " " , "." , grp_order2)

###########################################################################################
## 去除基因列，放到行名
rownames(lab_magic_RNA_meanbycluster) <- lab_magic_RNA_meanbycluster$gene
lab_magic_RNA_meanbycluster <- lab_magic_RNA_meanbycluster[,-1]

###########################################################################################
#####GSVA
tmp <- list()
for (i in unique(cell_research_dynamics$cell_type)){
  i <- list(subset(cell_research_dynamics,cell_type==i)$gene)
  tmp <- append(tmp,i)
  
}

names(tmp) <- unique(cell_research_dynamics$cell_type)
gsva_score <- gsva(as.matrix(lab_magic_RNA_meanbycluster),tmp,method="ssgsea",verbose=T,parallel.sz=40)

out_file <- paste0( out_path , "/" , out_name , "_gsva_score.csv" )
write.csv(gsva_score, out_file)

###########################################################################################
#pheatmap(R_studio)
cell_research_dynamics_gsva_score <- gsva_score
colnames(cell_research_dynamics_gsva_score) <- gsub( "Round.ElongateS" , "Round&ElongateS" , colnames(cell_research_dynamics_gsva_score) )
colnames(cell_research_dynamics_gsva_score) <- gsub( "Differenting.Differented" , "Differenting&Differented" , colnames(cell_research_dynamics_gsva_score) )
cell_research_dynamics_gsva_score <- cell_research_dynamics_gsva_score[,grp_order2]

p <- pheatmap(cell_research_dynamics_gsva_score,scale = "row",show_rownames = T ,
         #cutree_cols=7,cutree_rows = 7,
         cluster_rows = F,cluster_cols = F,clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
         cellwidth = 20, cellheight = 10)

out_file <- paste0( out_path , "/" , out_name , "_gsva_score.pdf" )
pdf(out_file)
print(p)
dev.off()